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Abstract 
We  describe  new  results  on  iterative  substructuring  methods. 
These  are  domain  decomposition  methods  for  which  the  subdomains,  also 
known  as  substructures,  do  not  overlap.  The  interaction  between  the 
subdomains,  to  satisfy  appropriate  continuity  requirements,  is  handled 
by  a  conjugate  gradient  method.  We  first  consider  the  case  of  two 
substructures  and  show  that  the  optimality  of  a  method  of  this  kind  can 
be  established  by  using  an  extension  theorem  for  finite  element  spaces. 
A  general  extension  theorem  is  then  formulated  for  general  conforming 
finite  elements  under  a  mild  restriction  on  the  triangulation.  The 
proof  of  this  result  extends  a  technique  introduced  by  Astrakhantsev. 
It  is  then  established  that  a  mechanism  for  global  transportation  of 
information  is  necessary  in  order  to  obtain  fast  convergence  in  the 
case  of  many  substructures.  In  the  case  of  a  simple  second  order 
problem  on  regions  in  the  plane,  an  alternative  derivation  is  given  of 
a  fast  mehod  introduced  by  Bramble,  Pasciak  and  Schatz.  Estimates  of 
the  rates  of  convergence  of  this  and  an  alternative  algorithm  are  also 
given. 
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Iterative  Substructuring  Methods;   The  General  Elliptic  Case 

by 
Olof  B.  Widlund 

1.  Introduction:  In  this  paper,  we  will  consider  a  family  of  domain 
decomposition  methods,  known  as  iterative  substructuring  methods,  which 
are  being  developed  to  solve  the  often  very  large  linear  systems  of 
algebraic  equations  that  arise  in  elliptic  finite  element  work.  These 
methods  differ  from  Schwarz's  alternating  algorithm  in  that  there  is  no 
overlap  of  the  subregions,  also  called  substructures,  into  which  the 
region  is  divided.  In  the  basic  form  of  these  algorithms,  the  finite 
element  or  finite  difference  approximations  of  the  elliptic  problem 
restricted  to  the  substructures,  are  solved  repeatedly  by  a  direct 
method  such  as  Gaussian  elimination  or,  in  special  cases,  a  fast 
Poisson  solver.  The  interaction  between  the  substructures,  to  satisfy 
various  continuity  requirements,  is  handled  by  an  iterative  method 
typically  a  preconditioned  conjugate  gradient  method.  For  a  discussion 
of  the  use  of  inexact  solvers  for  the  subproblems,  see  Bramble,  Pasciak 
and  Schatz  [7]  and  Dryja,  Proskurowski  and  Widlund,  [20,21]. 

The  idea  of  substructuring  of  large  discrete  elliptic  problems 
goes  back  at  least  to  the  early  1960's;  see  Przemienicki  [31].  In 
industrial  practice,  as  reflected  in  several  large  finite  element 
systems  such  as  NORSAM,  see  Bell,  Hatlestad,  Hausteen  and  Araldsen  [2], 
the  matrices  which  represent  the  interaction  between  the  substructures 
are  computed  and  they  are  then  factored  into  triangular  factors.  In 
such  programs  the  substructuring  idea  is  also  often  applied  recursively 
and  it  has  proven  a  very  useful  tool  for  organizing  software  systems 
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for  which  no  a  priori  upper  bound  on  the  problem  size  is  contemplated. 
We  note  that  in  such  contexts  it  is  important  that  the  modeling, 
factorization,  etc.,  of  substructures  can  be  saved  and  used  again 
together  with  one  or  a  few  remodeled  substructures  if  the  structure  is 
damaged.  Similarly  the  effects  of  a  redesign  of  part  of  one  of  these 
often  very  large  engineering  structures  can  be  modeled  at  a  much  more 
modest  expense. 

Since  the  arithmetic  work  of  algorithms  of  this  kind  almost 
invariably  grows  faster  than  linearly  in  the  number  of  degrees  of 
freedom,  a  divide  and  conquer  strategy,  implicit  in  the  entire 
substructuring  philosophy,  is  quite  attractive  provided  that  the 
interaction  phase  of  the  computation  can  be  carried  out  with  a  fast 
iterative  method.  We  note  that  Petter  Bj«$rstad  and  Anders  Hvidsten  of 
the  University  of  Bergen,  have  begun  an  ambitious  study  of  a  number  of 
large  industrial  problems  using  industrial  software  modified  to  allow 
the  use  of  iterative  substructuring  methods. 

Iterative  substructuring  methods  also  show  a  clear  promise  for 
parallel  computing.  We  can  expect  that  finite  element  codes  designed 
as  indicated  above  will  perform  quite  well  without  major  modifications 
on  parallel  systems  with  a  modest  number  of  processors.  It  can  also  be 
plausibly  argued  that  when  the  number  of  processors  grows,  bottlenecks 
might  occur  primarily  in  the  part  of  the  computation  that  involves  the 
interaction  between  the  substructures.  The  costs  of  this  phase,  both 
in  terms  of  arithmetic  and  communication,  can  be  considerably  decreased 
by  using  iterative  substructuring  methods.  No  firm  conclusions  can 
however   be   drawn  until   more  experience   has   been  gained  in  actual 


experiments.  We  note  that  Keyes  and  Gropp  [2^4]  have  already  undertaken 
such  a  study  on  an  Intel  Hypercube. 

The  case  of  two  or  a  few  substructures,  obtained  without  using 
intersecting  cuts,  is  by  now  quite  well  understood.  The  first  results 
on  optimal  methods  are  now  at  least  five  years  old;  see  Bjidrstad  and 
Widlund  [3,^4],  Bramble,  Pasciak  and  Schatz  [7],  Dryja  [16],  Lebedev 
[25],  Marchuk,  Kuznetscv  and  Matsokin  [26]  and  Widlund  [35,36].  A 
crucial  role  in  this  theory  is  played  by  an  extension  theorem  for 
finite  element  functions  which  until  recently  appears  to  have  been 
known  only  for  Lagrangian  elements  in  the  plane  and  the  H^-norm.  In 
section  2,  we  formulate  a  general  extension  theorem  for  general 
conforming  finite  element  spaces  in  an  arbitrary  number  of  dimensions. 
Our  proof,  given  in  Widlund  [36],  is  inspired  by  the  pioneering  work  of 
Astrakhantsev  [1]. 

In  section  3,  we  establish,  by  using  elementary  arguments,  that 
the  performance  of  iterative  substructuring  methods  must  decline  in  the 
absence  of  a  mechanism  for  global  transportation  of  information.  In 
that  section,  we  also  consider  the  simplest  case  of  problems  with 
intersecting  cuts,  namely  second  order  problems  in  planar  regions 
approximated  by  piecewise  linear  elements.  By  introducing  an  idea  from 
multigrid  methods,  we  derive  an  algorithm  which  turns  out  to  be 
identical  to  one  recently  introduced  by  Bramble,  Pasciak  and  Schatz 
[8].  An  alternative  algorithm,  see  Dryja  [17]  and  Dryja,  Proskurowski 
and  Widlund  [20,21],  is  also  described  and  estimates  are  given  of  the 
rate  of  convergence  of  these  two  algorithms. 

We  conclude  this  introduction  with  a  few  additional  comments  on 
the  literature.   Numerical  and  theoretical  work  on  the  case  of  many 
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substructures  without  intersecting  cuts  are  given  in  Dryja  and 
Proskurowski  [18,19].  Bramble,  Pasciak  and  Schatz  have  designed  and 
analyzed  methods  for  three  dimensional  problems,  see  [9].  Together 
with  Richard  Ewing  they  have  also  developed  related  methods  for  local 
mesh  refinements,  see  Bramble,  Ewing,  Pasciak  and  Schatz  [5]. 

There  has  also  been  an  active  interest  in  algorithms  of  this  kind 
at  Yale.  This  work  has  resulted  in  several  papers;  see  Chan  [10],  Chan 
and  Reasco  [11,12]  as  well  as  Keyes  and  Gropp  [2^4]  mentioned  earlier. 
Among  the  methods  considered  in  these  studies  is  one  due  to  Golub  and 
Mayers  [23].  For  a  discussion  of  a  number  of  different  methods, 
including  one  tested  very  early  by  Concus,  Golub  and  O'Leary  [14]  and 
several  suggested  by  Dihn,  Glowinski  and  Periaux  [15];  see  section  6  of 
Bj25rstad  and  Widlund  [4]. 

As  shown  for  example  in  Marchuk,  Kuznetsov  and  Matsdkin  [26]  and 
Widlund  [36],  these  methods  are  closely  related  to  capacitance  matrix 
methods,  known  in  the  Soviet  literature  as  fictitious  component  and 
fictitious  domain  algorithms.  The  resulting  systems  of  equations  have 
much  in  common  with  classical  potential  theory  developed  around  1900  by 
Fredholm  and  others;  see  Fredholm  [22],  Proskurowski  and  Widlund 
[29,30].  The  Soviet  literature,  with  which  this  author  unfortunately 
is  not  fully  acquainted,  also  traces  the  roots  of  continuous  analog  of 
the  substructuring  mehods  back  to  Poincar'e  and  Steklov,  also  active 
around  the  turn  of  the  century;  see  Lebedev  [25].  It  is  interesting  to 
note  that  Lebedev's  new  book  refers  to  composition  rather  than 
decomposition  methods.  In  a  case  of  large  engineering  structures  such 
as  oil  platforms  or  tankers,  this  is  indeed  a  more  natural  name  in  view 
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of  the  fact  that,  in  the  most  concrete  sense,  these  structures  must  be 
assembled  from  substructures. 
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2.  An  extension  theorem  for  finite  element  spaces.  We  begin  this 
section  by  considering  the  structure  of  the  finite  element  problems  and 
by  outlining  an  iterative  substructuring  algorithm.  We  do  this  in 
order  to  provide  motivation  for  the  development  of  an  extension 
theorem. 

Let  Q  be  a  Lipschitz  region  which  is  the  union  of  two  adjacent  but 
nonoverlapping  subregions  Q^  and  22  and  r3.  r3  is  the  intersection  of 
the  closures  of  fi^  and  ^z-  The  boundary  of  n  is  the  union  of  To  and  r^ 
on  which  essential  (Dirichlet)  and  natural  (Neumann)  boundary 
conditions  are  imposed.  The  boundary  of  Q  also  equals  T-^  U  Yz,  where 
the  boundary  of  n^  is  p^  u  r3,  etc.  A  linear,  selfadjoint  elliptic 
problem  of  order  22,  is  given  in  variational  form, 

ajj(u,v)  =  f(v)  ,   V  V  e  Vo(f2) 

where  Vq(j^)  is  the  subspace  of  the  Sobolev  space  H^(fi)  with  zero 
Dirichlet  data  on  Fq.  f(v)  is  a  linear  functional  obtained  from  the 
right  hand  side  of  the  elliptic  equation  and  the  natural  boundary 
conditions  by  using  a  Green's  formula.  We  assume  that  the  bilinear 
form  satisfies  the  following  standard  conditions, 

^n(u.v)  =  an(v,u) 
c|i"II^Q  i   an(u,u) 

|3fi(u,v)|  <  cIIuIIvqUvIIvq  , 
where  c  and  C  are  positive  constants.  We  can  then  use  the  Lax-Milgram 
theorem  to  prove  the  existence  of  a  unique  solution.  We  also  note  that 
the  discrete  elliptic  problem,  obtained  by  a  Galerkin  procedure  using  a 
conforming  finite  element  space,  will  have  a  positive  definite, 
symmetric  stiffness  matrix.   These  properties  follow  directly  from  the 


conditions  above  and  the  fact  that  the  finite  element  space  V*^  is  a 
subspace  of  V. 

Associated  with  the  finite  dimensional  space  V^  is  a  proper 
triangulation  of  the  subregions  f^^  '^2  ^^"'^  ^^  with  r^  following 
edges  of  triangles.  The  triangulation  can  be  quite  general  but  we 
insist  on  a  uniform  bound  on  h^/p\^,  where  h^  is  the  diameter  and  p^  the 
radius  of  the  largest  inscribed  sphere  of  the  element  K.  No  essential 
complications  arise  if  some  of  the  triangles  have  curved  edges  as  in 
isoparametric  methods,  cf.  Ciarlet  [13]. 

Following  Ciarlet  [13].  we  define  the  finite  element  space  in 
terms  of  a  pair  of  dual  bases  {(t>^}  and  {yj  }  satisfying  pi((f)j)  =  6ij  . 
Following  Strang  [33].  we  assume  that  the  basis  function  <^^  is 
supported  on  the  union  of  the  elements  associated  with  \i^  and  that  they 
all  are  uniform  of  degree  %  in  the  sense  that  there  exist  constants  Cg, 
independent  of  i  and  the  mesh  size,  such  that 

The  exponent  d^  is  often  the  degree  of  a  derivative  associated  with  the 
degree  of  freedom  y^.  it  can  be  defined  naturally  by  a  scaling 
argument  even  in  other  cases,  e.g.  when  \i^  represents  an  average  of  the 
function  or  one  of  its  derivatives  over  an  element  or  an  edge  of  an 

element.   The  natural  interpolation  problem  for  the  reconstruction  of 

,  di 

an  element  u^  ^  v  involves  the  quantities  h^    v^{u^) . 

The  linear  system  of  equations  corresponding  to  the  discrete 
elliptic  problem  has  the  form, 


Kx 


K^^  0  K13 
0  ^22  K23 
KJ3  k53  K33 


"w 
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■ 

XI 

r 

X2 

= 

b2 

''3 

b3 

(2.1) 


The  matrix  K11  represents  the  couplings  between  pairs  of  degrees  of 
freedom  associated  with  ^i  as  well  as  with  the  intersection  of  ffj  and 
f^l  ,  K13  the  couplings  between  pairs  belonging  to  this  set  and  those 
associated  with  r3  etc.  The  zero  blocks  in  the  (1.2)  and  (2.1) 
positions  reflect  the  fact  that  there  is  no  coupling  between  fii  and  Q2 
when  the  canonical  basis  functions  are  used.  The  elements  of  the 
stiffness  matrix  K  are  simply  obtained  as  k^^.  =  afj(  ^i  ,({)j  )  and  the  right 
hand  sides  from  f(<{>j)  and  the  Dirichlet  data. 

When  studying  iterative  substructuring  algorithms,  we  always 
assume  that  it  is  acceptable  to  solve  discrete  problem  on  the 
subregions,  with  some  appropriate  boundary  conditions  added  on  r3.  We 
can  therefore  concentrate  on  the  case  when  b-]  and  b2  are  zero,  since  we 
can  reduce  the  system  (2.1)  to  such  a  form  in  a  preliminay  step.  By 
block  Gaussian  elimination,  this  problem  is  reduced  to 


SX3  =  (K33  -  kI3KT]ki3  -  Ki3K2lK23)x3  =  b3 


(2.2) 


The  so  called  Schur  complement  is  given  by  S  =  S^^^  +  S^^' ,  where 
^  =  ^^2  ~  ^Tb^TI'^IS  is  the  Schur  complement  associated  with  a 
problem  on  the  subregion  U^    , 


'  Kn 

Kl3^ 

r       ] 

'  0' 

,  KT3 

vW 

.  ""3 

['I 

(2.3) 


Equation  (2.3)  is  the  finite  element  approximation  of  the  elliptic 
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problem  restricted  to  fi-,  ,  with  natural  boundary  conditions  added  on  r3. 
The  elements  of  K^^)  are  of  the  form  a^{,p^,(p^),  where  4)^  and  4,j  are 
basis  functions  associated  with  degrees  of  freedom  on  To.  We  note  that 
K33  in  (2.1)  naturally  can  be  written  as  K^3^  +  K^^^,  the  sum  of 
matrices  which  represent  the  contributions  to  the  integral  aj^(  (fjj^  ,(j)i ) 
from  f^i  and  ^2    >  respectively. 

It  has  been  shown,  see  e.g.  Bj«5rstad  and  Widlund  [3,'^]  or  Widlund 
[35],  that  the  condition  number  of  S  goes  to  infinity  when  the  problem 
size  increases.  A  preconditioner  should  therefore  be  found.  It  is  not 
surprising  that  S^^^  is  a  good  preconditioner  for  S  in  a  conjugate 
gradient  algorithm  to  solve  (2.2).  The  rate  of  convergence  is 
determined  by  the  stationary  values  of  the  generalized  Rayleigh 
quotient 

In  particular,  it  is  important  to  establish  that  these  eigenvalues  lie 
in  a  fixed  interval.  The  Rayleigh  quotient  is  bounded  from  below  by 
one  since  S  =  S^ ^ ^  +  S^^)  and  S^^^  is  positive  definite,  just  as  S^^^. 
An  upper  bound  of  (1+C)  can  be  established  by  showing  that 

x5s(2)x^  <  0x53^^^X3  ,   V  X3  {2.H) 

Since  b-i  and  b2  are  assumed  to  be  zero  in  equation  (2.1),  it  is 
natural  to  call  the  corresponding  finite  element  function  a  piecewise, 
discrete  harmonic  function.  It  is  easy  to  show  that  the  strain  energy 
x^Kx  of  such  a  function  is  equal  to  x^SXn  .  Similarly  xls^^^x^  and 
x^s'-'^^'x^  are  the  strain  energies  of  discrete  harmonic  functions  defined 
on  Qi  and  ^2  respectively,  with  common  boundary  values  on  r3, 
represented  by  Xo   . 
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Of  all  elements  in  h''(Q)  which  take  on  some  specific  values  on  dQ, 
the  solution  of  Laplace's  equation  provides  the  minimal  strain  energy. 
This  can  easily  be  shown  by  using  a  Green's  formula.  The  same  result 
holds  for  any  conforming  finite  elements.  Therefore,  to  establish  the 
bound  (2.4),  we  can  use  an  analog  of  a  well  known  extension  theorem  for 
Sobolev  spaces;  cf.  e.g.  Stein  [32].  We  assume  that  the  complement  of 
the  region  Q.  also  has  been  triangulated  in  an  equally  benign  way  and 
that  vh(n)  has  been  extended  to  V^CRn)  c  H^-CR'''). 

Theorem  1  .   Let  n  be   a  Lipschitz  region  and  let  v'^(fi)   be  a 

conforming  finite  element  space  spanned  by  basis  functions  satisfying 

the  uniformity  condition  given  above.   Then  there  exists  a  constant 

C(fi),  which  depends  on  the  Lipschitz  constant  of  the  boundary  of  fi,S,,n, 

the  choice  of  the  particular  finite  element  space  and  max(h  /n  )  but 

k  k 

not  on  the  dimension  of  V*^  and  a  linear  operator  E*^  which  maps  v'^(fl) 
into  v'^CR")  such  that, 

E%|^  =  Uh  ,      V  u^e   Vh(fi) 
and 

A  proof  is  given  in  VJidlund  [36].  We  note  that  its  key 
ingredients  are  the  extension  theorem  for  Sobolev  spaces  and  the  tools 
developed  by  Strang  [33]  to  approximate  general  elements  of  H^CR"^)  by 
finite  element  functions. 

The  proofs  given  previously  in  BjeJrstad  and  Widlund  [4]  and 
Bramble,  Pasciak  and  Schatz  [7]  were  limited  to  Lagrangian  and  linear 
finite  elements  respectively  to  plane  regions  and  to  the  case  K.  =  1  . 
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These  proofs,   virtually  identical,   use   Bramble-Hilbert ' s   lemma, 

Sobolev's  lemma  and  a  regularity  result  for  elliptic  problems  on 

regions  with  corners. 

We  note  that  it  is  easy  to  adopt  Theorem  1  to  obtain  a  bound  for 

the  iterative  substructuring  algorithm  considered  above.  We  need  to 

extend  a  finite  element  function  u  defined  on  Q^      and  vanishing  on 

h  ' 

Tq  n  Y]  into  a  function  on  fi  which  vanishes  on  rp.  First  extend  u^  by 
zero  into  bounded  subregions  of  the  complement  of  fi-j  which  are  bordered 
by  the  components  of  Fq.  This  step  leaves  the  strain  energy  unchanged 
and  we  can  now  use  the  theorem  given  above  to  extend  the  finite  element 
function  to  all  of  R".  Clearly  the  restriction  of  this  function  to  ^2 
provides  the  extension  necesary  co  establish  inequality  (2.4)  and  the 
optimality  of  the  conjugate  gradient  algorithm. 
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3.  Many  Substructures.  The  study  of  iterative  substructuring  methods 
for  problems  with  many  substructures  is  quite  important  since  we  wish 
to  find  out  if  these  algorithms  are  of  real  promise  for  parallel 
computers  with  many  processors.  We  are  of  course  primarily  interested 
in  algorithms  with  a  performance  which  does  not  deteriorate  very  much 
when  the  region  is  partitioned  into  an  increasing  number  of  subregions. 
The  algorithm  described  in  the  previous  section  works  very  well  for  two 
or  a  few  substructures  and  it  can  also  be  used  for  regions  partitioned 
into  any  number  of  strips.  However  in  each  iteration,  information  is 
exchanged  only  between  neighboring  substructures  across  the  interfaces. 
As  demonstrated  in  numerical  experiments  by  Dryja  and  Proskurowski 
[18],  [19],  the  rate  of  convergence  will  suffer  when  the  region  is 
divided  into  an  increasing  number  of  strips.  It  is  in  fact  easy  to  see 
that  a  residual  of  modestly  low  frequency  at  one  end  of  a  region  will 
affect  the  error  everywhere  and  that  therefore  the  number  of  steps 
until  convergence  of  any  iterative  method  must  in  general  exceed  the 
number  of  strips  if  the  interaction  is  solely  through  next  neighbors. 
A  similar  argument  is  valid  for  general  regions.  Speedup  can  only  be 
obtained  by  introducing  a  mechanism  for  global  transportation  of 
information  across  the  region. 

A  solution  is  provided  by  considering  an  idea  from  multigrid 
theory;  see  Yserentant  [37].  Let  us  for  simplicity  consider  only 
piecewise  linear  elements  and  Poisson's  equation,  in  which  case 


aj^(u,v)  =  /  Vu'Vv  dx  =  (u,v)  ,    , 


and  a  region  in  the  plane  divided  into  triangular  substructurs  Q^    , 


each  of  which  in  turn  is  divided  into  triangular  elements.  We  only 
consider  the  Dirichlet  case  and  use  the  standard  notation  for  the  norms 
and  seminorms  of  Sobolev  spaces.  As  in  section  2,  we  avoid  very  thin 
triangles  in  both  of  these  triangulations.  The  typical  diameters 
associated  with  the  two  partitions  are  H  and  h  respectively. 

As  explained  in  section  2,  we  can  work  throughout  with  piecewise, 
discrete  harmonic  functions.  Any  such  function  Uj^  can  be  partitioned 
into 

"h  =  iH^h  +  (uh  -  iH^^h)  (3.1) 

where  Ij^  is  the  interpolation  operator  associated  with  the  coarse  mesh, 
cf.  Yserentant  [37].  It  is  easy  to  show  that  the  piecewise  linear 
function  l^u^  as  well  as  u^-lH^^h  ^""^  piecewise,  discrete  harmonic. 
Since  the  second  term  in  (3.1)  vanishes  at  all  vertices  of  the  coarse 
mesh,  it  represents  higher  frequencies;  the  pitch  of  a  drum  increases 
substantially  if  held  fixed  at  a  number  of  points.  A  preliminary 
preconditioner  is  introduced  in  terms  of  the  bilinear  form, 

(lHUh,lHVh)j^1(jj)  +  (i^h-lH%.Vh-lHVh)^l(j^)  •  (3.2) 

The  resulting  linear  system  of  equations  is  block  diagonal  if  we  work 
with  the  standard  basis  functions  on  the  coarse  mesh  and  the  standard 
ones  for  the  fine  mesh  associated  with  all  the  vertices  except  those 
which  are  also  a  vertex  of  a  large  triangle.  One  of  the  blocks  thus 
corresponds  exactly  to  the  original  problem  solved  on  the  coarse  mesh. 
If  this  system  of  equations  is  solved  directly,  then  it  provides  the 
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necessary,  quite  modest  amount  of  global  information.  The  second  term 
should  be  modified  since  by  Itself  it  gives  rise  to  a  linear  system 
which  is  almost  as  difficult  as  the  original  one. 

Before  this  is  done  it  is  important  to  note  that  for  v^^  =  u^  the 
expression  (3.2)  is  bounded  from  below  by  J.  luhl^i  .  This  follows 
from  the  triangle  inequality.   The  first  term  lln^ihl^i    is  bounded  by 


'H^(n) 


"^1 '"hi  1,00^   s"*^  ^^   finally  obtain 


\^h^h\l^,,,  i   C(i  *  iog(H/h))  luhl^i   .      (3.3) 


HMn)         °      '  '"HUn) 


The  same  type  of  bound  can  also   easily   be   established   for 

l^^h-lHUhl  1,  V   by  using  the  triangle  inequality.  The  estimate  (3.3) 
"   'H' (fl) 

follows  from  an  inequality  for  finite  element  functions  in  the  plane 
that  does  not  have  a  direct  continuous  counterpart, 


see  Bramble  [5],  Bramble,  Pasciak  and  Schatz  [8],  Rukhovets  and 
Oganesyan  [27],  Thomee  [3^]  or  Yserentant  [37].  This  result  is  often 
given  only  for  piecewise  linear  functions  but  it  also  holds  for  more 
general  finite  elements.  Note  that  log(H/h)  measures  the  number  of 
refinement  steps  needed  to  go  from  the  substructure  level  to  the  fine 
triangulation. 

The  estimate  (3- ■4)  cannot  be  improved,  but  the  Lp-term  can  be 
removed  for  the  case  at  hand  by  noting  that 

Ih(uh+const . )  I  1,    =  I  luUh+const .  I  ■,  ,  =  llHUhl  i, 

"  "       'hUq^)   '  "  "^       'hU^^)   '  "  "'H4n.) 

and  that  therefore  Poincare's  inequality  can  be  applied. 
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We  can  therefore  conclude  that  the  preconditioner  given  by  (3.2) 
gives  a  conjugate  gradient  mehod  with  a  condition  number  bounded  by 
C(l  +  log(H/h))  and  that  it  therefore  would  converge  in 
C  log(l/e)  ( 1 +log(H/h) ) '' '^^  steps.  If  we  destroyed  the  mechanism  for 
global  transportation  of  information  by  replacing  the  coarse  finite 
element  model,  defined  by  the  first  term  of  (3.2),  by  a  diagonal 
matrix,  then  the  condition  number  estimate  would  increase  by  a  factor 
(1/H)2  and  the  number  of  iteration  steps  to 
C  log(l/£)  (1+log(H/h))l/2(i/H). 

As  in  the  previous  section,  we  can  associate  a  strain  energy  with 

each  individual  substructure.   The  estimates  given  below  will  be  for 

individual   substructures  and  the  final  results  will  simply  be  obtained 

by  adding  the  contributions  from  all  of  them.   In  Q^   the  finite  element 

function  Wj^  =  uh-lH^ih  depends  only  on  its  boundary  values  since  it  is 

discrete  harmonic  in  n^.  we  will  obtain  our  final  preconditioner  by 

replacing   Iwhl^i,     ,  one  of  the  contributions  from  Qi      to  the 

preliminary  preconditioner  (3.2),  by  the  square  of  a  norm  of  these 

boundary  values.   Before  we  proceed,  we  introduce  a  second  optional 

preconditioner  for  the  case  of  two  substructures.   It  has  been 

established   in  BjeJrstad  and  Widlund  [^4]  that  S^^^,  defined  in 

section  2,  is  spectrally  equivalent  to  J.   This  operator  is  defined  in 

terms  of   its  square  J^  which  is  a  finite  element  discretization  of 

-(A)2  on  To   with  Dirichlet  conditions  at  the  end  points.   In  this 
ds       -> 

construction,  we  use  the  mesh  and  finite  element  space  obtained  by 
restricting  the  given  finite  element  space  to  the  curve  r^.  As  shown 
in  BjeJrstad  and  Widlund  [^],  the  quadratic  form  x^Jx^  is  equivalent  to 


whe 
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I    ,?  ,   ,2         ^3  r  l^h(>^(s))|2    |uh(x(s))|2 

r-e   X,3  is  the  length  of  r3.   The  seminorm  l^hl  i/2(p   i^  given  by 

^3  ^3  |u(x(s))-u(x(s'))|^  ,3  ,3.  ^         (3.5,) 
'"'^lHl/2(r3)  =  (5  f)  |s-s'|2 

Let  3Q^  =  ^  Fij  ,  where  Ti  j  =  jii  n  ^j  i  '^  J  .  ^n  edge  of  the 
triangular  substructure  Q^.  Since  we  know  how  to  work  with  the 
operator  J  and  Wj^  vanishes  at  the  vertices  of  the  substructures,  it  is 
natural  to  replace  |wi.|2      in  the  preconditioner  by 


This  is 


The  expression  (3-6)  provides  an  upper  bound  for  l^hluW 
so  since  it  is  known  that  any  element  in  H^^^Cfy)  can  be  extended  by 
zero  to  the  rest  of  dn^  with  the  resulting  function  in  H  On^).  The 
HUn^)-norm  of  the  discrete  harmonic  function  with  these  boundary 
values  can  be  bounded  by  the  H^^2(r.,)  norm  of  its  trace  on  r^j.  This 
follows  from  the  extension  theorem  given  in  Section  2;  for  details  see 
Bj(Z5rstad  and  Widlund  [^].  Since  w^  can  be  written  as  the  sum  of  three 
functions  of  this  kind  we  obtain, 


|2,  ,    <  C  I  IIwhI|2wo,_    .  (3.7) 


^'<H.o'']   l'"^"Hi^^(r,,: 


To  obtain  an  upper   bound  for   Z  M^hll  i/2c     '   ^®  ^^^  ^^^ 
definition  (3.5).   We  use  the  trace  theorem  and  find  that 
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Since  the  seminorm  of  the  left  hand  side  does  not  change  if  we  add  a 
constant  to  w^,  we  can  use  Poincare's  inequality  and  the  bound  obtained 
previously  for  |w|,^|       to  obtain, 

It  is  also  straightforward  to  estimate  the  second  term  of   (3.5a). 
By  using  that  w^  vanishes  at  the  vertices  of  ^i,  we  find  that 

«-ij  |wj^(x(s))|2     hij  |^,^(jj(3))|2     i^.    |„^(j^(3))|2 

/  ds  =  /  ds  +  /  ds 

OS  6       s         h^j      s 


|wv^(x(h;  j))|2  hij 


i     ^  d^   -    "''hi  1^00(5^.)  1°S  (dij/hij)(3.9) 
i   C(l  -  log(H/h))l|wj^||2^^^^  <  c(l  .  log(H/h))  Muhll^„(^^  • 

Here  h^j  is  the  distance  to  the  first  mesh  point  along  T^j   and  l^j    the 
length  of   r^j.    The   last   inequality   follows   easily   since 

We  obtain  from  O.'^),  (3.8)  and  (3.9)  that 

J  """""J^^lrij)  -  <:('*iog(H/h))2(|u,|2,^_^^^  ,  ,/„2  ll-hil^2,,^,)  .(3.10) 

and  finally,   since  w^  =  (i-ij^)uh  =  (  I-Ih)  (uh+const . )  ,   the  following 
bounds  are  established  by  another  application  of  Poincare's  inequality. 


=l"hl^,,„.,  S  |lH"hlJ,(,.,  *I  l|uh-lH"«ll^,^2,r^,, 


C(l  +  log(H/h)2)  |u^,|2 

"'hUq^) 
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Using  standard  technical  tools,  we  have  just  provided  is  a 
simplified  proof  of  the  main  result  in  Bramble,  Pasciak  and  Schatz  [8]. 

Theorem  2.  The  preconditioned  conjugate  gradient  method  which 
uses  the  preconditioner  defined  by 


(lHUh.lHVh),i(,)  -  ci  .1.  ((uh-lHUh.  Vh-lHVh))^j/2(, 


j) 


converges  in  C  log(1/£)*(1  +  log(H/h))  steps.   Here  ci  is  any  constant 
on  the  order  of  one. 

The  same  tools  can  be  used  to  prove  an  equally  strong  bound  for  a 
natural  extension  of  the  method  discussed  in  section  2.  This  algorithm 
is  described  more  fully  in  Dryja  [17]  and  Dryja,  Prokurowski  and 
Widlund  [20], [21].  In  brief,  we  divide  the  substructures  into  two 
disjoint  sets.  We  denote  the  corresponding  unions  by  il^''  and  Q^^' . 
These  sets  will  play  the  same  role  as  Q-^  and  fip  i"  ^^^  algorithm 
described  in  section  2.  We  must  assume  that  the  substructuring  is  such 
that  if  two  substructures  Q^  and  fij  have  a  common  edge  r^j  =  fi^  n  fij  . 
then  fi^  and  Qj  do  not  belong  to  the  same  subset.  The  set  Q]  '  is  thus 
a  union  of  the  closure  of  certain  substructures,  and  these 
substructures  intersect  only  at  vertices  of  substructures.  We  can 
formulate  a  Neumann  problem  on  n^^^  and  it  is  possible  to  reduce  the 
corresponding  large  linear  system  of  equations  to  linear  systems 
corresponding  to  individual  substructures  by  using  Lagrange 
multipliers;  for  details  see  Dryja  [17]  and  Dryja,  Proskurowski  and 
Widlund  [20],  [21].  Just  as  described  in  section  2,  the  solution  of 
the  linear  system  which  corresponds  to  the  preconditioner  amounts  to 
solving  a  Neumann  problem  on  n^^^  followed  by  a  Dirichlet  problem  on 
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fi(2)  using  the  boundary  values  obtained  from  the  Neumann  problem  as 
Dirichlet  data.  In  the  case  at  hand,  Q^^'  is  the  disjoint  union  of  a 
number  of  substructures.  These  problems  are  therefore  independent  and 
can  be  solved  in  parallel.  We  also  note  that  global  transportation  of 
information  is  assured  by  the  fact  that  Q^"^'    is  connected. 

As  in  section  2,  the  rate  of  convergence  is  determined  by  an  upper 
bound  for 


x^Sx^/x^S^^^x^  , 


i.e.  by  bounding  the  strain  energy  associated  with  ^^^'  in  terms  of 
the  energy  of  Q^^'.  This  can  be  accomplished  by  estimating  the  strain 
energy  of  one  substructure  belonging  to  the  second  set  by  the  strain 
energy  of  those  neighbors  which  belong  to  the  first  set. 

In  our  proof  we  again  partition  up^  as  in  (3.1).  The  strain  energy 
attributable  to  Ij^uj^  can  be  handled  straightforwardly  by  using 
inequality  (3-3)  for  individual  subregions.  From  inequality  (3-10) 
follows  that 

M-h-lHUhM2,l/2(,..^^C(l  .log(H/h)2)  |u,l2^^^^  . 

i.e.  the  square  of  the  Hj,^^-norm  of  the  boundary  data  of  u^-lp^uj^  on 
one  edge  of  a  region  can  be  bounded  by  the  strain  energy  of  one  of  its 
next  neighbors  at  the  expense  of  a  factor  (1  +  log(H/h)2).  The  use  of 
estimate  (3.7)  now  completes  the  proof  of  the  following  theorem. 

Theorem  3.    The   preconditioned   conjugate   gradient   method 
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introduced  by   Dryja,   Proskurowski   and   Widlund   converges   in 
C  log(l/e)(l  +  log(H/h))  steps. 

We  note  in  conclusion  that  this  theory  has  been  extended  to  more 
general  families  of  finite  element  spaces;  see  Dryja  [17]  and  Dryja, 
Proskurowski  and  Widlund  [21 ]. 
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